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Jacobi-type iterative algorithms for the eigenvalue decomposition, singular value decomposition, and Takagi 
factorization of complex matrices are presented. They are implemented as compact Fortran 77 subroutines in a 
freely available library. 



1. Introduction 

This note describes a set of routines for the 
eigenvalue decomposition, singular value decom- 
position, and Takagi factorization of a complex 
matrix. Unlike many other implementations, the 
current ones are all based on the Jacobi algo- 
rithm, which makes the code very compact but 
suitable only for small to medium-sized problems. 

Although distributed as a library, the routines 
are self-contained and can easily be taken out of 
the library and included in own code, removing 
yet another installation prerequisite. Owing to 
the small size of the routines (each about 3 kBytes 
source code) it is possible, in fact quite straight- 
forward, to adapt the diagonalization routine to 
one's own conventions rather than vice versa. 

2. Mathematical Background 

2.1. Eigenvalue Decomposition 

The eigenvalue decomposition of a nonsingular 
matrix A £ (C" x ™ takes the form 

UAU- 1 =diag(<7i,...,<7„), Oi e C. (1) 

The eigenvalues Oi and transformation matrix U 
can be further characterized if A possesses certain 
properties: 

• A = A^ (Hermitian): U~ l = W, cr t S H, 

• A = A T (symmetric): U^ 1 = U T . 

2.2. Singular Value Decomposition 

The singular value decomposition (SVD) can 
be applied to an arbitrary matrix A £ ([™x« : 

V*AW* =diag(<7i,...,<7 fi ), (2) 



n = min(m, n) , Cj ^ . 

V consists of orthonormal row vectors, i.e. is also 
unitary for m = n. 

2.3. Takagi Factorization 

The Takagi factorization [1,2] is a less known 
diagonalization method for complex symmetric 
matrices A = A T £ <C nx, \ 



XJ*AXJ^ =diag(<7i,...,cr„), 
t/- 1 = £/•+, Oi^O. 



(3) 



Although outwardly similar to the eigenvalue de- 
composition of a Hermitian matrix, it is really the 
special case of an SVD with V — W* , as it applies 
even to singular matrices. Note also that the left 
and right factors, U* and IP, are in general not 
inverses of each other. 

One might think that the Takagi factorization 
is merely a scaled SVD. For example, the matrix 
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has the SVD 

V T diag(o-i,cr 2 )VF 



(4) 
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which can indeed be scaled to yield 
J7 T diag(<7i, <7 2 )U = 



(6) 
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But consider the matrix 
A = 



which has the SVD 
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whereas its Takagi factorization is 
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Although occurring less frequently than the 
eigenvalue decomposition, the Takagi factoriza- 
tion does have real applications in physics, e.g. in 
the diagonalization of mass matrices of Majorana 
fermions. 

3. Jacobi Algorithm 

The Jacobi algorithm [3] consists of iteratively 
applying a basic 2x2 diagonalization formula un- 
til the entire n x n matrix is diagonal. It works 
in several 'sweeps' until convergence is achieved. 
In each sweep it rotates away the non-zero off- 
diagonal elements using the 2x2 algorithm. Ev- 
ery such rotation of course creates other non-zero 
off-diagonal elements. It can be shown, however, 
that the sum of the absolute values of the off- 
diagonal elements is reduced in each sweep. More 
precisely, the Jacobi method has quadratic con- 
vergence [4]. 

Convergence is in most cases achieved in 5—10 
sweeps, which for an n x n matrix translates into 
(10— 20)n 3 multiply-add operations to obtain the 
eigenvalues only and (15— 30)n 3 operations in- 
cluding the eigenvectors ([5], cf. also Sect. 10). 
This compares with |n 3 + 30n 2 operations for the 
Householder /QL algorithm when just the eigen- 
values are sought and |n 3 + 3n 3 when also the 
eigenvectors are needed. 

For large n, the Jacobi algorithm is thus not 
the most efficient method. Nevertheless, for small 



to medium-sized problems the Jacobi method is 
a strong competitor, in particular as it has the 
following advantages: 

• It is conceptually very simple and thus very 
compact. 

• It delivers the eigenvectors at little extra cost. 

• The diagonal values are accurate to machine 
precision and, in cases where this is mathemat- 
ically meaningful, the vectors of the transfor- 
mation matrix are always orthogonal, almost 
to machine precision. 

For the various diagonalization problems dis- 
cussed before, only the core 2x2 diagonalization 
formula changes, whereas the surrounding Jacobi 
algorithm stays essentially the same. 

The famous Linear Algebra Handbook gives an 
explicit implementation of the Jacobi algorithm 
for real symmetric matrices [4] , taking particular 
care to minimize roundoff errors through mathe- 
matically equivalent variants of the rotation for- 
mulas. The present routines are closely patterned 
on this procedure. For the Takagi factorization, 
the use of the Jacobi algorithm was first advo- 
cated in two conference papers [6,7] which give 
only few details, however. The two-sided Jacobi 
version of the SVD is used in [8] . Ref . [9] gives a 
more thorough overview of literature on the Ja- 
cobi method. 

4. The 2x2 Formulas 

4.1. Eigenvalue decomposition 

Using the ansatz 



U 



ci iici 



(10) 



y -t 2 c 2 c 2 
the equation UA = diag(<7i, a 2 ) U becomes 

cri = A n +t x A 2 i = A 22 + — A12 , (11) 

11 

0-2 = A u - — A 21 = A 22 - t 2 A 12 . (12) 
t 2 



Solving for t\ and t 2 yields 

A12 ^ A 2 i 

t l — H — 1 l 2 — 



A + D 



A + D' 



(13) 
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A = -(An - A 22 ) , 
D = ±yj A 2 + A 12 A 21 . 



(14) 
(15) 



For the numerical stability it is best to choose the 
sign of D which gives ii j2 the larger denominator. 
This corresponds to taking the smaller rotation 
angle (< 7r/4). The diagonal values are 



a 1 =A 11 +5, _ A 12 A 21 
<r 2 = A 22 -5, ~ A + D 



(16) 



In order that U smoothly becomes unitary as A 
becomes Hcrmitian, we choose 



C l = C 2 = : , 

VI + *1*2 

which guarantees a unit determinant. 

4.2. Takagi Factorization 

Substituting the unitary ansatz 



(17) 



U 



c = 



c tce lip 
-tce-' lv c 

' , teJR, 



VTTf 1 

into U*A = diag(<7i, a 2 ) U and introducing 



(T x = e'^cri 



An = e'^A 



n , 



(18) 
(19) 

(20) 



o 2 = c -^a 2 , A 22 = c-^A 22 , (21) 
we arrive at 

ffi = in + *Ai2 = A 22 + jA 12 , (22) 

a 2 = In - jA 12 = A 22 - tA 12 . (23) 

Comparing with Eqs. (11) and (12), the solution 
can be read off easily: 



t 



112 



A + D 
A = \{A xl -A 22 ) 



d = ±Ja* + a*. 



(24) 
(25) 
(26) 



Again it is best for numerical stability to choose 
the sign of D which gives the larger denominator 
for t. The diagonal values become 



a 1 =A 11 +tA 12 e^, 
g 2 = A 22 -tA 12 e iv . 



(27) 
(28) 



The assumption t £ 1R fixes the phase ip. It re- 
quires that Ai 2 and A have the same phase, i.e. 
A = (real number) • A\ 2 . Since both e 1(p and its 
conjugate appear in A, we try the ansatz 

= aA 12 + (3A* 12 (29) 

and choose coefficients to make the A\ 2 term in 



A oc (aAn - f3*A 22 )Ai 2 + 
(/?A n - a*A 22 )A{ 2 



(30) 



vanish. This is achieved by a = A n and /3 = A 22 
which also makes the coefficient of Ai 2 real. Thus, 



A* U A 12 + A 22 A\ 2 



\A* U A 12 + A 22 A* 12 \ 
5. Singular Value Decomposition 



(31) 



We insert unitary parameterizations for the left 
and right transformation matrices X = V, W, 



X = 



cx = 



cx 



tx c x 



-t* x c x c x 



(32) 
(33) 



into V* A = diag(cri, a 2 ) W and by eliminating 
<7i ; 2 arrive at 

A 12 + A 22 ty = (An + A 21 ty)t w , (34) 
An + A 22 t* w = (An + A 12 t* w )t v . (35) 

The solutions are evidently related through ex- 
change of the off-diagonal elements. Explicitly, 



M v 



tY ~ A v + D V ' 
M v = A* n A 21 + A 22 A* 12 , 

Ay = ^(|An| 2 -|A 22 | 2 + 

|Ai 2 | 2 -L4 2 i| 2 ), 



(36) 
(37) 

(38) 
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D v = ±yjA*, + \M v \*, (39) 

tw = W\ Al2 ^ A21 ■ (40) 

Dy and Dyy share the same sign, which is chosen 
to yield the larger set of denominators for better 
numerical stability. 

The singular values become 

<ri = —(A 11 +fyA 21 ), (41) 

Cw 

o 2 = —{A 22 -t v A 12 ), (42) 

Cw 

If A e (D mx ™ is not a square matrix, we con- 
sider two cases: 

For m > n, we make A square by padding it 
with zero-columns at the right. For zero right col- 
umn, A12 = A 22 = 0, Eq. (42) guarantees that it 
is d 2 that vanishes. That is, the above Jacobi ro- 
tation never moves a singular value into the zero- 
extended part of the matrix. All singular values 
automatically end up as the first min(m, n) diag- 
onal values of the Jacobi-rotated A. 

For m < n, we apply this algorithm to A T and 
at the end exchange V and W. This is the least 
involved solution, as A has to be copied to tem- 
porary storage for zero-extension anyway. 

6. Installation 

The Diag package can be downloaded from the 
URL http : //www . f eynarts . de/diag. After un- 
packing the tar file, the library is built with 

. /configure 
make 

To compile also the Mathematica executable, one 
needs to issue "make all" instead of just "make." 
The generated files are installed into a platform- 
dependent directory with "make install" and at 
the end one can do a "make clean" to remove 
intermediate files. 

The routines in the Diag library allocate space 
for intermediate results according to a preproces- 
sor variable MAXDIM, defined in diag.h. This ef- 
fectively limits the size of the input and output 
matrices but is necessary because Fortran 77 of- 
fers no dynamic memory allocation. Since the Ja- 
cobi algorithm is not particularly suited for large 



problems anyway, the default value of 16 should 
be sufficient for most purposes. 

7. Description of the Fortran Routines 

The general convention is that each matrix is 
followed by its leading dimension in the argument 
list, i.e. the m in A(m,n). In this way it is pos- 
sible to diagonalizc submatrices with just a dif- 
ferent invocation. Needless to add, the leading 
dimension must be at least as large as the corre- 
sponding matrix dimension. 

7.1. Hermitian Eigenvalue Decomposition 

Hermitian matrices are diagonalized with 

subroutine HEigensystem(n, A.ldA, 

d, U.ldU, sort) 
integer n, IdA, ldU, sort 
double complex A(ldA,n), U(ldU,n) 
double precision d(n) 

The arguments are as follows: 

• n (input), the matrix dimension. 

• A (input), the matrix to be diagonalized. Only 
the upper triangle of A needs to be filled and it 
is further assumed that the diagonal elements 
are real. Attention: the contents of A are not 
preserved. 

• d (output), the eigenvalues. 

• U (output), the transformation matrix. 

• sort (input), a flag that determines sorting of 
the eigenvalues: 

= do not sort, 

1 = sort into ascending order, 
— 1 = sort into descending order. 

The 'natural' (unsortcd) order is determined by 
the choice of the smaller rotation angle in each 
Jacobi rotation. 

7.2. Symmetric Eigenvalue Decomposition 

The second special case is that of a complex 
symmetric matrix: 
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subroutine SEigensystem(n, A,ldA, 

d, U,ldU, sort) 
integer n, IdA, ldU, sort 
double complex A(ldA,n), U(ldU,n) 
double complex d(n) 

The arguments have the same meaning as for 
HEigensystem, except that A's diagonal elements 
are not assumed real and sorting occurs with re- 
spect to the real part only. 

7.3. General Eigenvalue Decomposition 

The general case of the eigenvalue decomposi- 
tion is implemented in 

subroutine CEigensystem(n, A, IdA, 

d, U.ldU, sort) 
integer n, IdA, ldU, sort 
double complex A(ldA,n), U(ldU,n) 
double complex d(n) 

The arguments are as before, except that A has 
to be filled completely. 

7.4. Takagi Factorization 

The Takagi factorization is invoked in almost 
the same way as SEigensystem: 

subroutine TakagiFactor (n, A, IdA, 

d, U.ldU, sort) 
integer n, IdA, ldU, sort 
double complex A(ldA,n), U(ldU,n) 
double precision d(n) 

The arguments are as for SEigensystem. Also 
here only the upper triangle of A needs to be filled. 

7.5. Singular Value Decomposition 

The SVD routine has the form 

subroutine SVD(m, n, A, IdA, 

d, V.ldV, W.ldW, sort) 
integer m, n, IdA, ldV, ldW, sort 
double complex A(ldA,n) 
double complex V(ldV,m), W(ldW,n) 
double precision d(min(m,n)) 

with the arguments 

• m, n (input), the dimensions of A. 

• A (input), the m x n matrix of which the SVD 
is sought. 



• d (output), the singular values. 

• V (output), the min(m, n) xm left transformation 
matrix. 

• W (output), the min(m,n) x n right transforma- 
tion matrix. 

• sort (input), the sorting flag with values as 
above. 

8. Description of the C Routines 

The C version consists merely of an include 
file CDiag.h which sets up the correct interfacing 
code for using the Fortran routines. In particu- 
lar the usual problem of transposition 1 between 
Fortran and C is taken care of. 

The arguments arc otherwise as for Fortran. To 
treat complex numbers uniformly in C and C++, 
CDiag.h introduces the new double_complex 
type which is equivalent to complex<double> in 
C++ and to struct { double re, im; } in C. 

C's syntax unfortunately does not allow the 
declaration of variable-size matrices as function 
arguments, thus it is not possible for CDiag.h to 
set up the prototypes to directly accept C matri- 
ces without compiler warnings. To simplify mat- 
ters, CDiag.h defines the macro Matrix which is 
used as in 

double_complex A [5] [3] , V [3] [5] , W[3] [3] ; 
double d[3] ; 

SVD(5, 3, Matrix(A) , d, 

Matrix (V), Matrix (W), 0); 

This is equivalent to using an explicit cast and 
passing the leading dimension, i.e. 

SVD (5, 3, (double_complex *)A,3, d, 
(double_complex *)V,5, 
(double_complex *)W,3, 0); 

Compilation should be done using the included 
fee script, i.e. by replacing the C compiler with 
fee, for example 

1 C uses row- major storage for matrices, whereas Fortran 
uses column-major storage, i.e. a matrix is a vector of row 
vectors in C, and a vector of column vectors in Fortran. 
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fee -Iinclude myprogram . c -Llib -ldiag 

The fee script is configured and installed dur- 
ing the build process and automatically adds the 
necessary flags for linking with Fortran code. 

9. The Mathematica Interface 

The Mathematica version may not seem as use- 
ful as the Fortran library since Mathematica al- 
ready has perfectly functional eigen- and singular 
value decompositions. The Takagi factorization 
is not available in Mathematica, however, and 
moreover the interface is ideal for trying out, in- 
teractively using, and testing the Diag routines. 

The Mathematica executable is loaded with 

Install ["Diag"] 

and makes the following functions available: 

• HEigensystem [A] computes the eigenvalue de- 
composition {d , U} of the Hermitian matrix A 
such that U.A == DiagonalMatrix [d] .U. 

• SEigensystem [A] computes the eigenvalue de- 
composition {d , U} of the symmetric matrix A 
such that U.A == DiagonalMatrix [d] .U. 

• CEigensystem [A] computes the eigenvalue de- 
composition {d , U} of the general matrix A such 
that U.A == DiagonalMatrix [d] .U. 

• TakagiFactor [A] computes the Takagi factor- 
ization {d, U} of the symmetric matrix A such 
that 

Conjugate [U] .A == DiagonalMatrix [d] .U. 

• SVD [A] computes the singular value decom- 
position {V , d , W} of the matrix A such that 
Conjugate [V] .A == DiagonalMatrix [d] .W. 

Note that these routines do not check whether 
the given matrix fulfills the requirements, e.g. 
whether it is indeed Hermitian. 

10. Timings 

The following plot shows the time for diago- 
nalizing a random matrix of various dimensions. 
Note that the abscissa is divided in units of the di- 
mension cubed; this accounts for the anticipated 



scaling behaviour of the Jacobi algorithm, hence 
the curves appear essentially linear. 

The absolute time values should be taken for 
orientation only, as they necessarily reflect the 
CPU speed. For reference, the numbers used in 
the figure above were obtained on an AMD64-X2 
CPU running at 3 GHz. Each point is actually 
the average from diagonalizing 10 6 random ma- 
trices, to reduce quantization effects in the time 
measurement. 




The average number of sweeps needed to diag- 
onalize the 10 6 random matrices to machine pre- 
cision is plotted in the next figure. 



# sweeps 




11. Summary 

The Diag package contains Fortran subroutines 
for the eigenvalue decomposition, singular value 
decomposition, and Takagi factorization of a com- 
plex matrix. The Fortran library is supplemented 
by interfacing code to access the routines from 
C/C++ and Mathematica. 

The routines are based on the Jacobi algorithm. 
They are self-contained and quite compact, thus 
it should be straightforward to use them outside 
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of the library. All routines are licensed under the 
LGPL. 
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